*Title: Delivering remote learning using a low-tech solution: Evidence from an RCT during school closures  
*Author: Liang Choon Wang, Michael Vlassopoulos, Asad Islam, Hashibul Hassan 

*These codes generate figure 4-6 included in the main text

***************************************************************************************
*Stata preparation
***************************************************************************************
*net install grc1leg2, from("http://digital.cgdev.org/doc/stata/MO/Misc") replace 
*ssc install schemepack, replace
*ssc install coefplot
graph set window fontface "Times New Roman"
set scheme white_tableau
set scheme plotplain

***************************************************************************************
**Figure 4. Treatment effects on the standardized learning outcomes 
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

*Temporary Excel file
global control_list 		child_age baseline_literacy_score baseline_numeracy_score child_pvt_tuition father_edu_years mother_edu_years hh_mem family_income television smartphone hs_land reli_dummy


putexcel set std_coeff_gender.xlsx, replace 
putexcel A1 = "regressor" B1 = "beta" C1 = "se" D1 = "DoF" E1 = "treat_arms" F1 = "gender"

putexcel 	A2 	= "Total score" 	A3 	= "Literacy" 		A4 	= "Numeracy" 
putexcel 	A5 	= "Total score" 	A6	= "Literacy" 		A7 	= "Numeracy"
putexcel 	A8 	= "Total score" 	A9 	= "Literacy" 		A10 = "Numeracy"
putexcel 	A11 = "Total score" 	A12	= "Literacy" 		A13 = "Numeracy"
putexcel 	A14 = "Total score" 	A15	= "Literacy" 		A16 = "Numeracy"
putexcel 	A17 = "Total score" 	A18	= "Literacy" 		A19 = "Numeracy"

local std_outcome_vars endline_total_score_std endline_literacy_score_std endline_numeracy_score_std 

local j = 2

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list child_gen i.child_grade_std if endline_completed == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[standard_treat]  C`j' = _se[standard_treat] D`j' = (e(df_r)) E`j' = "Standard" F`j' = "Pooled"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list child_gen i.child_grade_std if endline_completed == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[extended_treat]  C`j' = _se[extended_treat] D`j' = (e(df_r)) E`j' = "Extended" F`j' = "Pooled"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1 & child_gen == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[standard_treat]  C`j' = _se[standard_treat] D`j' = (e(df_r)) E`j' = "Standard" F`j' = "Boy"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1 & child_gen == 0, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[standard_treat]  C`j' = _se[standard_treat] D`j' = (e(df_r)) E`j' = "Standard" F`j' = "Girl"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1 & child_gen == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[extended_treat]  C`j' = _se[extended_treat] D`j' = (e(df_r)) E`j' = "Extended" F`j' = "Boy"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1 & child_gen == 0, r cluster(VILLAGE_ID)
putexcel set std_coeff_gender.xlsx, modify
putexcel B`j' = _b[extended_treat]  C`j' = _se[extended_treat] D`j' = (e(df_r)) E`j' = "Extended" F`j' = "Girl"
local j = `j' + 1
}

import excel "std_coeff_gender.xlsx", firstrow clear 
erase std_coeff_gender.xlsx

gen 		levelN = .
replace 	levelN = 1	if treat_arms == "Standard"
replace 	levelN = 2	if treat_arms == "Extended"

gen			type = 1 if regressor == "Total score" & gender == "Pooled"
replace	type = 2 if regressor == "Literacy" & gender == "Pooled"
replace	type = 3 if regressor == "Numeracy" & gender == "Pooled"
replace	type = 4 if regressor == "Total score" & gender == "Boy"
replace	type = 5 if regressor == "Literacy" & gender == "Boy"
replace	type = 6 if regressor == "Numeracy" & gender == "Boy"
replace	type = 7 if regressor == "Total score" & gender == "Girl"
replace	type = 8 if regressor == "Literacy" & gender == "Girl"
replace	type = 9 if regressor == "Numeracy" & gender == "Girl"

gen 		regOrder = .
replace 	regOrder = 1		if gender == "Pooled"
replace 	regOrder = 2 	if gender == "Boy"
replace 	regOrder = 3 	if gender == "Girl"

sort type regOrder levelN
gen order = _n

foreach h in 5 9{
gen max9`h'=. // generating variable for upper bound of CI (95% & 99%)
gen min9`h'=. // generating variable for lower bound of CI (95% & 99%)
}

foreach h in max min{

	replace `h'99 = invttail(89, 0.005)
	replace `h'95 = invttail(89, 0.025)

	}

foreach h in 5 9{
gen max9`h'b = .
gen min9`h'b = .
}

* upper bound
foreach h in 5 9{
	replace max9`h'b = beta + se*max9`h'
}

*lower bound
foreach h in 5 9{
	replace min9`h'b = beta - se*min9`h'
}

order regressor min*95*b max*95*b min*99*b max*99*b beta se order

sort order 

keep beta m*b levelN regressor order gender treat_arms

order beta m*b levelN regressor order

sort order

forval j=1/2{
mkmat beta-levelN if levelN==`j', matrix(A`j')
matrix A`j'2=A`j''
mat list A`j'2
 }

coefplot matrix(A12[1,]) matrix(A22[1,]), ci((2 3) (4 5))  msymbol(t) msize(medsmall) ///
xline(0, lwidth(medium) lcolor(gray) lpattern(solid)) ///
xscale(range(-0.30 1)) ///
ylabel(, labsize(vsmall) nogrid) ///
xlabel(-0.30(0.1)1, labsize(vsmall)) /// 
xtitle("Effect size in standard deviations of the control group", size(vsmall)) ///
coeflabels( 	r1 = "Total score" ///
						r2 = "Literacy" ///
						r3 = "Numeracy" ///
						r4 = "Total score" ///
						r5 = "Literacy" ///
						r6 = "Numeracy" ///
						r7 = "Total score" ///
						r8 = "Literacy" ///
						r9 = "Numeracy") ///
headings(r1 = "{bf: Panel A: Pooled Sample}" r4 = "{bf: Panel B:  Boys Sample}" r7 = "{bf: Panel C: Girls Sample}", ///
gap(0.5) labsize(vsmall)) ciopts(recast(rcap rspike) lwidth(medium thin)) ///
legend(order(3 "Standard" 6 "Extended") region(col(none)) size(vsmall)  rows(2) position(10) ring(0)) ///
graphreg(color(gs16)) aspectratio(1.2)  mlabel(string(beta,"%9.2f")) mlabposition(2)  mlabsize(tiny) name(Figure_4, replace)

graph export "Figure 4.emf", replace
*graph export "Figure 4.eps", replace

****************************************************************************************
*Figure 5. Treatment effects on secondary outcomes
****************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

*Temporary Excel file
global control_list 	child_age child_gen baseline_literacy_score baseline_numeracy_score child_pvt_tuition father_edu_years mother_edu_years hh_mem family_income television smartphone hs_land reli_dummy

putexcel set std_coeff_others.xlsx, replace 
putexcel A1 = "regressor" B1 = "beta" C1 = "se" D1 = "DoF" E1 = "treat_arms"

putexcel 	A2 	= "Overall impulsivity" 				A3 = "Schoolwork impulsivity" 		A4 	= "Interpersonal impulsivity" ///
			A5 	= "Grit score" 							A6 = "Growth mindset" 				A7	= "Affective empathy (Contagion)" ///
			A8 	= "Cognitive empathy (Understanding)" 										A9 	= "Prosocial motivation (Support)" ///
			A10 = "Leadership" 							A11 = "Communication" 				A12	= "Planning" ///
			A13	= "Emotional symptoms" 					A14 = "Conduct problem" 			A15 = "Hyperactivity" ///
			A16	= "Peer problem" 						A17 = "Prosocial"

putexcel 	A18 = "Overall impulsivity" 				A19 = "Schoolwork impulsivity" 	A20	= "Interpersonal impulsivity" ///
			A21	= "Grit score" 							A22 = "Growth mindset" 				A23	= "Affective empathy (Contagion)" ///
			A24 = "Cognitive empathy (Understanding)" 	A25 = "Prosocial motivation (Support)" ///
			A26 = "Leadership" 							A27 = "Communication" 				A28	= "Planning" ///
			A29	= "Emotional symptoms" 					A30 = "Conduct problem" 			A31 = "Hyperactivity" ///
			A32	= "Peer problem" 						A33 = "Prosocial"

local std_outcome_vars 	impulsivity_overall_std impulsivity_schoolwork_std impulsivity_interpersonal_std grit_std gms_std /// 
						affective_empathy_std cognitive_empathy_std prosocial_empathy_std leadership_std ///
						communication_std planning_std sdq_emotion_std sdq_conduct_std sdq_hyper_std ///
						sdq_peer_std sdq_prosocial_std

local j = 2

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_others.xlsx, modify
putexcel B`j' = _b[standard_treat]  C`j' = _se[standard_treat] D`j' = (e(df_r)) E`j' = "Standard"
local j = `j' + 1
}

foreach v of local std_outcome_vars {
regress `v' standard_treat extended_treat $control_list i.child_grade_std if endline_completed == 1, r cluster(VILLAGE_ID)
putexcel set std_coeff_others.xlsx, modify
putexcel B`j' = _b[extended_treat]  C`j' = _se[extended_treat] D`j' = (e(df_r)) E`j' = "Extended"
local j = `j' + 1
}

import excel "std_coeff_others.xlsx", firstrow clear 

erase "std_coeff_others.xlsx"

gen 	levelN = 1 if treat_arms == "Standard"
replace levelN = 2 if treat_arms == "Extended"

gen type=.
replace type = 1 	if inlist(regressor, "Leadership")
replace type = 2	if inlist(regressor, "Communication")
replace type = 3	if inlist(regressor, "Planning")
replace type = 4 	if inlist(regressor, "Overall impulsivity")
replace type = 5 	if inlist(regressor, "Schoolwork impulsivity")
replace type = 6 	if inlist(regressor, "Interpersonal impulsivity")
replace type = 7 	if inlist(regressor, "Grit score")
replace type = 8 	if inlist(regressor, "Growth mindset")
replace type = 9 	if inlist(regressor, "Affective empathy (Contagion)")
replace type = 10 	if inlist(regressor, "Cognitive empathy (Understanding)")
replace type = 11 	if inlist(regressor, "Prosocial motivation (Support)")
replace type = 12 	if inlist(regressor, "Emotional symptoms")
replace type = 13 	if inlist(regressor, "Conduct problem")
replace type = 14 	if inlist(regressor, "Hyperactivity")
replace type = 15 	if inlist(regressor, "Peer problem")
replace type = 16 	if inlist(regressor, "Prosocial")

gen 		regOrder = 1 	if regressor == "Leadership"
replace 	regOrder = 2	if regressor == "Communication"
replace 	regOrder = 3	if regressor == "Planning"
replace 	regOrder = 4 	if regressor == "Overall impulsivity"
replace 	regOrder = 5 	if regressor == "Schoolwork impulsivity"
replace 	regOrder = 6 	if regressor == "Interpersonal impulsivity"
replace 	regOrder = 7 	if regressor == "Grit score"
replace 	regOrder = 8 	if regressor == "Growth mindset"
replace 	regOrder = 9 	if regressor == "Affective empathy (Contagion)"
replace 	regOrder = 10 	if regressor == "Cognitive empathy (Understanding)"
replace 	regOrder = 11 	if regressor == "Prosocial motivation (Support)"
replace 	regOrder = 12 	if regressor == "Emotional symptoms"
replace 	regOrder = 13 	if regressor == "Conduct problem"
replace 	regOrder = 14 	if regressor == "Hyperactivity"
replace 	regOrder = 15 	if regressor == "Peer problem"
replace 	regOrder = 16 	if regressor == "Prosocial"

sort type regOrder levelN
gen order = _n

foreach h in 5 9{
gen max9`h'=. // generating variable for upper bound of CI (95% & 99%)
gen min9`h'=. // generating variable for lower bound of CI (95% & 99%)
}

foreach h in max min{

	replace `h'99 = invttail(89, 0.005)
	replace `h'95 = invttail(89, 0.025)

	}

foreach h in 5 9{
gen max9`h'b = .
gen min9`h'b = .
}

* upper bound
foreach h in 5 9{
	replace max9`h'b = beta + se*max9`h'
}

*lower bound
foreach h in 5 9{
	replace min9`h'b = beta - se*min9`h'
}

order regressor min*95*b max*95*b min*99*b max*99*b beta se order

sort order 

keep beta m*b levelN regressor order

order beta m*b levelN regressor order

sort order

forval j=1/2{
mkmat beta-levelN if levelN==`j', matrix(A`j')
matrix A`j'2=A`j''
mat list A`j'2
 }

coefplot matrix(A12[1,]) (matrix(A22[1,])), ci((2 3) (4 5)) msymbol(t) msize(medsmall) /// 
xline(0, lwidth(medium) lcolor(gray) lpattern(solid)) ///
xscale(range(-1 1)) ///
ylabel(, labsize(vsmall) nogrid) ///
xlabel(-1(0.25)1, labsize(vsmall)) xtitle("Effect size in standard deviations of the control group", size(vsmall)) ///
coeflabels( 	r1 = "Leadership" ///
						r2 = "Communication" ///
						r3 = "Planning" ///
						r4 = "Overall impulsivity" ///
						r5 = "Schoolwork impulsivity" ///
						r6 = "Interpersonal impulsivity" ///
						r7 = "Grit score" ///
						r8 = "Growth mindset" ///
						r9 = "Affective empathy (Contagion)" ///
						r10 = "Cognitive empathy (Understanding)" ///
						r11 = "Prosocial motivation (Support)" ///
						r12 = "Emotional symptoms" ///
						r13 = "Conduct problem" ///
						r14 = "Hyperactivity" ///
						r15 = "Peer problem" ///
						r16 = "Prosocial" ///
						) ///
headings(r1 = `""{bf: Panel A:  Leadership, communication,}" "{bf: & planning skills}""' ///
r4 = "{bf: Panel B: Noncognitive skills}" r12 = "{bf: Panel C:  Behavioural difficulties}", gap(0.5) labsize(vsmall)) ///
ciopts(recast(rcap rspike) lwidth(medium thin)) graphreg(color(gs16)) aspectratio(1.2) ///
legend(order(3 "Standard" 6 "Extended") region(col(none)) size(vsmall)  rows(2) position(4) ring(0)) ///
title("") mlabel(string(beta,"%9.2f")) mlabsize(tiny) mlabposition(2) mlabvposition(beta) mlabgap(0) name(Figure_5, replace)

graph export "Figure 5.emf", replace
*graph export "Figure 5.eps", replace

****************************************************************************************
*Figure 6. Heterogeneity in assessment test performance, by baseline numeracy & literacy score, family income, parents' education quartiles  
****************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

*****************************************************************************************************************
*Barplot - income_quartile 
*****************************************************************************************************************
drop if endline_completed == 0
xtile family_income_quartile = family_income, n(4)

collapse (mean) mean = endline_total_score (sd) sd = endline_total_score (count) n = endline_total_score, by(treatment_arms family_income_quartile)
generate hi = mean + invttail(n-1,0.025)*(sd / sqrt(n))
generate lo = mean - invttail(n-1,0.025)*(sd / sqrt(n))

generate income = treatment_arms    if family_income_quartile == 1
replace  income = treatment_arms+4  if family_income_quartile == 2
replace income = treatment_arms+8	if family_income_quartile == 3
replace  income = treatment_arms+12  if family_income_quartile == 4

sort income
list income family_income_quartile treatment_arms, sepby(family_income_quartile)

twoway 	(bar mean income if treatment_arms==1, fcolor(black) lcolor(black) lpattern(solid)) ///
		(bar mean income if treatment_arms==2, fcolor(none) lcolor(black) lpattern(solid)) ///
		(bar mean income if treatment_arms==3, fcolor(none) lcolor(black) lpattern(dash_dot)) ///
		(rcap hi lo income, lcolor(black)), ///
		legend(order(1 "Standard" 2 "Extended" 3 "Control" 4 "95% CI") rows(1) pos(6)) ///
		xlabel(2 "1{sup:st} quartile" 6 "2{sup:nd} quartile" 10 "3{sup:rd} quartile" 14 "4{sup:th} quartile", noticks) ///
		 xtitle("Family income") ytitle("Total score") name(income, replace)

****************************************************************************************
*graph - barplot - result 
****************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

drop if endline_completed==0

egen baseline_result = rowtotal(baseline_literacy_score baseline_numeracy_score)
xtile baseline_result_quartile = baseline_result, n(4)

collapse (mean) mean = endline_total_score (sd) sd = endline_total_score (count) n = endline_total_score, by(treatment_arms baseline_result_quartile)
generate hi = mean + invttail(n-1,0.025)*(sd / sqrt(n))
generate lo = mean - invttail(n-1,0.025)*(sd / sqrt(n))

generate result = treatment_arms    if baseline_result_quartile == 1
replace  result = treatment_arms+4  if baseline_result_quartile == 2
replace result = treatment_arms+8	if baseline_result_quartile == 3
replace  result = treatment_arms+12  if baseline_result_quartile == 4

sort result
list result baseline_result_quartile treatment_arms, sepby(baseline_result_quartile)

twoway	(bar mean result if treatment_arms==1, fcolor(black) lcolor(black) lpattern(solid)) ///
		(bar mean result if treatment_arms==2, fcolor(none) lcolor(black) lpattern(solid)) ///
		(bar mean result if treatment_arms==3, fcolor(none) lcolor(black) lpattern(dash_dot)) ///
		(rcap hi lo result, lcolor(black)), ///
		legend(order(1 "Standard" 2 "Extended" 3 "Control" 4 "95% CI") rows(1) pos(6)) ///
		xlabel(2 "1{sup:st} quartile" 6 "2{sup:nd} quartile" 10 "3{sup:rd} quartile" 14 "4{sup:th} quartile", noticks) ///
		xtitle("Baseline numeracy & literacy score") ytitle("Total score") name(result, replace)

*********************************************************************************
*graph - barplot - Parental Education
*********************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

drop if endline_completed==0

egen parentsedu = rowtotal(father_edu_years mother_edu_years)
xtile parentsedu_quartile = parentsedu, n(4)

collapse (mean) mean = endline_total_score (sd) sd = endline_total_score (count) n = endline_total_score, by(treatment_arms parentsedu_quartile)
generate hi = mean + invttail(n-1,0.025)*(sd / sqrt(n))
generate lo = mean - invttail(n-1,0.025)*(sd / sqrt(n))

generate	pedu = treatment_arms    	if parentsedu_quartile == 1
replace  	pedu = treatment_arms + 4  	if parentsedu_quartile == 2
replace 	pedu = treatment_arms + 8 	if parentsedu_quartile == 3
replace  	pedu = treatment_arms + 12 	if parentsedu_quartile == 4

sort pedu
list pedu parentsedu_quartile treatment_arms, sepby(parentsedu_quartile)

twoway 	(bar mean pedu if treatment_arms==1, fcolor(black) lcolor(black) lpattern(solid)) ///
		(bar mean pedu if treatment_arms==2, fcolor(none) lcolor(black) lpattern(solid)) ///
		(bar mean pedu if treatment_arms==3, fcolor(none) lcolor(black) lpattern(dash_dot)) ///
		(rcap hi lo pedu, lcolor(black)), ///
		legend(order(1 "Standard" 2 "Extended" 3 "Control" 4 "95% CI") rows(1) pos(6)) ///
		xlabel(2 "1{sup:st} quartile" 6 "2{sup:nd} quartile" 10 "3{sup:rd} quartile" 14 "4{sup:th} quartile", noticks) ///
		 xtitle("Parents' education in years") ytitle("Total score")  name(pedu, replace)

*combine
grc1leg2 result income pedu,  ycommon  iscale(*.9) cols(1) name(Figure_6, replace)
graph display, ysize(12) xsize(7)

graph export "Figure 6.emf", replace
*graph export "Figure 6.eps", replace

clear 